
##突变图

setwd("E:\\1.Ovarian_Treg\\1.picture\\12.MUT")

library(RColorBrewer)
library(ComplexHeatmap)
## 突变二值矩阵
mut <- read.table("mutMatrix.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
#转换为1-0矩阵
st=function(x){
  ifelse(x=="Mutation",1,0)
}
mut=apply(mut, 2, st)

## 突变频数
tmb <- read.table("TMB.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)

## 突变签名
mutsig <- read.table("rs_tcga.txt",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
mutsig$risk=ifelse(mutsig$riskScore>median(mutsig$riskScore),"high","low")
mutsig=mutsig[,c("riskScore","risk"),drop=F]
# 或者换成下面这行，读取FigureYa110mutationSignature的输出文件mutsig.weightMatrix.txt
# mutsig <- read.table("mutsig.weightMatrix.txt", row.names = 1, header = T) 


## 亚型数据
subt <- read.table("rs_tcga.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
subt$risk=ifelse(subt$riskScore>median(subt$riskScore),"high","low")
subt=subt[,c("risk"),drop=F]
#colnames(subt)="CMOIC"
# 设置亚型颜色
clust.col <- ggsci::pal_npg("nrc")(10)#,"#EC7D21"
blue   <- "#5bc0eb"
red    <- "#f25f5c"

# 处理突变签名数据
#mutsig <- mutsig[,c("Signature.1","Signature.2","Signature.5","Signature.13")] # 文章中用到3种类型的signature， SBS1 (age-related), SBS2 and SBS13 (APOBEC activity-related) and SBS5 (ERCC2 mutation-related)
mutsig$APOBEC <- mutsig$riskScore # APOBEC相关的签名由签名2和3叠加
#mutsig$Cluster <- subt[rownames(mutsig),"Cluster"] # 添加亚型结果
#mutsig$Cluster.NUM <- ifelse(mutsig$Cluster=="C1",1,
#ifelse(mutsig$Cluster=="C2",2,3)) 
mutsig <- mutsig[order(mutsig$risk,decreasing = T),] # 确定整个热图的排序，按照亚型升序以及APOBEC降序排列
# 挑选要展示的基因
mutgene <- read.table("mut_gene.txt")[,1]

# 制作oncoprint的输入数据
sam=intersect(rownames(mutsig),colnames(mut))
mut=mut[,sam]
mutsig=mutsig[sam,]
onco.input <- mut[mutgene,rownames(mutsig)]
onco.input[onco.input == "1"] <- "Mutated" # 二值矩阵中1记为突变
onco.input[onco.input != "Mutated"] <- "" # 非“突变”给予空值
alter_fun = list(
  background = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = "#dcddde", col = "#dcddde"))
  },
  Mutated = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = "#A60000", col = "#A60000")) 
  }
)
col = c("Mutated" ="#A60000") # 突变颜色，注意这里只给了主图像的图例

my_ann <- subt[rownames(mutsig),,drop = F]
#my_annotation2 = HeatmapAnnotation(df = my_ann[,1,drop=F])
my_annotation1 = HeatmapAnnotation(df = my_ann[,1,drop=F], 
                                   col = list(risk = c("high" = clust.col[1],
                                                       "low" = clust.col[2])))
# 突变主区域的上部注释（突变负荷柱状图）
tmb$TMB=log2(tmb$TMB+1)
colnames(tmb)="log10TMB"
top_anno <- anno_barplot(as.numeric(tmb[rownames(mutsig),"log10TMB"]),
                         border = FALSE,
                         gp = gpar(fill = "#3379B4",border =NA,lty="blank"), 
                         height = unit(2.5, "cm"))

# 突变主区域的上部注释（突变签名柱状图）
tmp <- mutsig[,c("riskScore"),drop=F] # 只取和APOBEC有关的签名
#tmp$Others <- 1 - rowSums(tmp) # 计算其他签名的比例
top_anno2 <- anno_barplot(as.numeric(log2(tmp$riskScore+1)),
                          border = FALSE,
                          gp = gpar(fill = "#8491B4FF", 
                                    border = NA, # 无边框
                                    lty = "blank"),
                          height = unit(2, "cm")) # 高度

tmp <- as.data.frame(t(mut[mutgene,rownames(mutsig)]))

mut.order <- names(sort(colSums(tmp),decreasing = T)) # 根据突变频数高低排序展示突变的顺序
tmp$risk <- subt[rownames(tmp),"risk"]
pct <- NULL # 计算各个基因突变的百分比
for (i in mut.order) {
  tmp1 <- tmp[,c(i,"risk")]
  tmp1 <- as.data.frame.array(table(tmp1[,1],tmp1$risk))[2,]/sum(tmp1[,1])
  pct <- rbind.data.frame(pct,tmp1)
}
rownames(pct) <- mut.order

# 添加右侧百分比堆积柱状图
right_anno <- anno_barplot(as.matrix(pct),
                           which = "row",
                           border = FALSE,
                           gp = gpar(fill = clust.col,border=NA,lty="blank"), 
                           bar_width = 0.6,
                           width = unit(1.8, "cm"),
                           height = unit(1, "cm"))
###底部注释

library(circlize)
clinical=read.table("TMEscores.txt", header=T, sep="\t", check.names=F, row.names=1)

clinical=clinical[colnames(onco.input),]
my_annotation2 <- HeatmapAnnotation(
  ESTIMATEScore=clinical$ESTIMATEScore,
  StromalScore=clinical$StromalScore,
  ImmuneScore=clinical$ImmuneScore,
  TumorPurity=clinical$TumorPurity)




p2 <- oncoPrint(onco.input[mut.order,rownames(my_ann)], # 排序的突变矩阵
                alter_fun = alter_fun,  # 主区域的函数，包括各单元格大小、背景颜色等等
                col = col, # 突变颜色
                bottom_annotation = NULL, # 无底部注释
                top_annotation = c(HeatmapAnnotation(TMB = top_anno), # 顶部第一个注释：TMB
                                   my_annotation1, # 顶部第二个注释：亚型
                                   my_annotation2),
                # bottom_annotation=column_ha,# 顶部第三个注释：突变签名
                column_order = rownames(my_ann), # 样本的排序，根据突变签名的顺序
                right_annotation = rowAnnotation(PCT = right_anno), # 右侧堆叠柱状图注释
                show_pct = T, # 展示左侧的百分比
                column_title = "", # 不显示主题
                show_heatmap_legend = T, # 展示图例
                column_split = my_ann$CMOIC, # 根据亚型切分热图
                column_title_gp = gpar(fontsize = 8),
                row_names_gp = gpar(fontsize = 8),
                column_names_gp = gpar(fontsize = 8))
draw(p2)
library(ggpubr)
library(reshape2)
library(limma)
library(ggplot2)
riskFile="rs_tcga.txt"      #风险文件
TMBFile="TMB.txt"             #肿瘤突变负荷文件
setwd("E:\\1.Ovarian_Treg\\1.picture\\12.MUT")      #设置工作目录

#读取输入文件
TMB=read.table(TMBFile, header=T, sep="\t", check.names=F, row.names=1)      #读取肿瘤突变负荷文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)    #读取风险文件
risk$risk=ifelse(risk$riskScore>median(risk$riskScore),"high","low")
#合并数据
TMB=as.matrix(TMB)
TMB[TMB>quantile(TMB,0.975)]=quantile(TMB,0.975)
sameSample=intersect(row.names(TMB), row.names(risk))
TMB=TMB[sameSample,,drop=F]
risk=risk[sameSample,,drop=F]
data=cbind(risk, TMB)
data=data[,c("riskScore", "risk", "TMB")]
data$TMB=log2(data$TMB+1)
#data$riskScore=log2(data$riskScore+1)
#设置比较组
data$risk=factor(data$risk, levels=c("low", "high"))
risk=levels(factor(data$risk))
comp=combn(risk, 2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
#绘制箱线图
boxplot=ggplot(data = data,aes(x = risk, #分组列名
                               y = TMB, #连续变量列名
                               fill = risk))+ #按分组填充颜色
  scale_fill_manual(values = c("#4DBBD599","#E64B3599")) + #用自定义颜色填充
  geom_violin(alpha = 0.4, position = position_dodge(width = .75), 
              size = 0.8, color="black") +
  geom_boxplot(notch = FALSE, outlier.size = -1, 
               color="black", lwd=0.8, alpha = 0.7) +
  geom_point(shape = 21, size=2, # 点的性状和大小
             position = position_jitterdodge(), # 让点散开
             color="black", alpha = 1) +
  theme_classic() + 
  ylab("log2(TMB+1)") +
  xlab("Risk") +
  theme(axis.text.x = element_text(hjust = 1, size = 12,face = "bold.italic"),
        #panel.border = element_rect(colour = "black", fill=NA, size=0.2),
        axis.ticks = element_line(size=0.2, color="black"),
        axis.ticks.length = unit(0.2, "cm"),
        legend.position = "none",
        axis.title = element_text(size = 15,face = "bold.italic"),
        axis.text = element_text(size = 10)) + 
  # 如果不要组间比较就注释掉下面这行
  stat_compare_means(comparisons = my_comparisons)
pdf(file="boxplot.pdf",width=5,height=4.5)
print(boxplot)
dev.off()

#相关性图形
data=data[order(data$riskScore,decreasing = T),]
data=data[6:nrow(data),]
p1=ggplot(data=data, aes(x=riskScore, y=TMB))+ 
  xlab("Risk score") + ylab("log2(TMB+1)")+
  geom_point(color = "#988d7b") + 
  geom_smooth(method = "lm", 
              formula = y ~ x,
              fill = "#273359", 
              color = "#00aeef", 
              alpha = 0.8) + 
  theme_bw()+
  theme(
    # 去除网格线：
    panel.grid = element_blank(),
    # 修改坐标轴标签
    axis.title = element_text(face = "bold.italic"),
    # 标题居中：
    plot.title = element_text(hjust = 0.5, size = 10)
  )+
  stat_cor(method = 'spearman', aes(x =riskScore, y =TMB))
#相关性图形
pdf(file="TMBcor.pdf", width=3, height=3)
print(p1)
dev.off()

##survival

#引用包
library(survival)
library(survminer)

riskFile="rs_tcga.txt"      #风险文件
tmbFile="TMB.txt"             #肿瘤突变负荷文件

#读取输入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)    #读取风险文件
risk$risk=ifelse(risk$riskScore>median(risk$riskScore),"high","low")
tmb=read.table(tmbFile, header=T, sep="\t", check.names=F, row.names=1)      #读取肿瘤突变负荷文件

#合并数据
sameSample=intersect(row.names(tmb), row.names(risk))
tmb=tmb[sameSample,,drop=F]
risk=risk[sameSample,,drop=F]
data=cbind(risk, tmb)

#获取肿瘤突变负荷最优的cutoff
res.cut=surv_cutpoint(data, time = "futime", event = "fustat", variables =c("TMB"))
cutoff=as.numeric(res.cut$cutpoint[1])
tmbType=ifelse(data[,"TMB"]<=cutoff, "L-TMB", "H-TMB")
scoreType=ifelse(data$risk=="low", "low risk", "high risk")
mergeType=paste0(tmbType, "+", scoreType)

#定义生存分析函数
bioSurvival=function(surData=null, outFile=null){
  diff=survdiff(Surv(futime, fustat) ~ group, data=surData)
  length=length(levels(factor(surData[,"group"])))
  pValue=1-pchisq(diff$chisq, df=length-1)
  if(pValue<0.001){
    pValue="p<0.001"
  }else{
    pValue=paste0("p=",sprintf("%.03f",pValue))
  }
  fit <- survfit(Surv(futime, fustat) ~ group, data = surData)
  #print(surv_median(fit))
  
  #绘制生存曲线
  bioCol=c( "#00A087FF" ,"#4DBBD5FF","#E64B35FF","#3C5488FF")
  bioCol=bioCol[1:length]
  surPlot=ggsurvplot(fit, 
                     data=surData,
                     conf.int=F,
                     pval=pValue,
                     pval.size=6,
                     legend.title="Risk-TMB",
                     legend.labs=levels(factor(surData[,"group"])),
                     font.legend=10,
                     legend = c(0.8, 0.8),
                     xlab="Time(years)",
                     break.time.by = 1,
                     palette = bioCol,
                     surv.median.line = "hv",
                     risk.table=F,
                     cumevents=F,
                     risk.table.height=.25,
                     ggtheme = theme(panel.border = element_rect(fill=NA,
                                                                 color="black", 
                                                                 size=1, 
                                                                 linetype="solid"))+
                       
                       theme(panel.grid = element_blank()),
                     title = "Survival curve", # 添加标题
                     font.main = c(16, "bold.italic", "black"), # 设置标题字体大小、格式和颜色
                     font.x = c(14, "bold.italic", "black"), # 设置x轴字体大小、格式和颜色
                     font.y = c(14, "bold.italic", "black"), # 设置y轴字体大小、格式和颜色
                     font.tickslab = c(12, "plain", "black"),# 设置坐标轴刻度字体大小、格式和颜色
                     linetype = "solid")
  
  pdf(file=outFile, width=5.5, height=4.8, onefile = FALSE)
  print(surPlot)
  dev.off()
}

#调用函数,绘制肿瘤突变负荷的生存曲线
data$group=tmbType
bioSurvival(surData=data, outFile="TMB.survival.pdf")

#调用函数,绘制肿瘤突变负荷联合病人风险的生存曲线
data$group=mergeType
bioSurvival(surData=data, outFile="TMB-risk.survival.pdf")
save(data,file = "TMB.group.rdata")

